-- require "geometry"
local nl=require "luaNumeric";
local unity = require "luaUnity";
require "luaFDTD3D"


--===========================================================================--
--初始化参数
local RadomeSolver={

	scan_angles = {0,10};
	beam_types = {'SUM','DIFF' }; --和波束，DIFF是差波束
	polars = {'HH','VV'};   --极化，有交叉极化，罩子的耦合比较小，-20dB左右
	polarHH = true;
	polarVV = false;
	beamSUM = true;
	beamDIFF = false;
	
	epsr = 3.0;
	tans = 0.001;
	tans_fit_freq = 10;
	
	posx_antenna = 350;  --天线位置，天线的绝对位置，假设天线放在x轴上
	radius_antenna = 100;  ---天线口面半径， 
	ant_min_freq = 8;
	ant_max_freq = 12;
	step_angle = 31;
	min_angle = 0;
	max_angle = 60;
	is_phasearray = false;
	
	farbox_min_freq = 10;
	farbox_max_freq = 11;
	farbox_num_freq = 3;
	mesh_freq = 12;
	accuracy = -30;
	max_steps = 10000;
	num_threads = 20;
	cal_air = 0;
	cal_nearfields = true;

	

	--theta_antenna = 0;           --天线与x轴夹角，天线口面与x轴夹角
	--polarization = "HH";            --极化，有交叉极化，罩子的耦合比较小，-20dB左右
	--theta_beam = 20;             --天线指向，天线自身的偏转

	set_settings_str = Gallop.set_settings_str;
	get_settings_str = Gallop.get_settings_str;
};
--===========================================================================--

--===========================================================================--
--创建表和创建元表
function RadomeSolver:new (o)
	o = o or {};
	self.__index = self;
	setmetatable(o,self)
	--self.update_parameters_from_ui();
	o:regist_funcs_to_gallop();
	o.init_strs = o:get_settings_str();
	return o;
end
--===========================================================================--


--===========================================================================--
--注册全局函数
function RadomeSolver:regist_funcs_to_gallop()
	Gallop.init_funcs['RadomeSolverApp'] = function()
		local str = self.init_strs;
		assert( type('str')=='string', 'RadomeSolverApp Plugin init string error' );
		if str ~= '' then
			Gallop.RadomeSolverApp:set_settings_str(str)
		end
		Gallop.RadomeSolverApp:update_parameters_to_ui();
	end;
	
	Gallop.load_funcs['RadomeSolverApp'] = function()
		local str = Geometry.get_setting( 'RadomeSolverApp.Settings' )
		if str ~= '' then
			Gallop.RadomeSolverApp:set_settings_str(str)
		end
		Gallop.RadomeSolverApp:update_parameters_to_ui();
	end;
	
	Gallop.presave_funcs['RadomeSolverApp'] = function()
		Gallop.RadomeSolverApp:update_parameters_from_ui()
		local str = Gallop.RadomeSolverApp:get_settings_str();
		Geometry.set_setting( 'RadomeSolverApp.Settings', str );
	end;
end
--===========================================================================--

--===========================================================================--
--GUI交互函数
function RadomeSolver:update_parameters_to_ui()
	--print( 'update parameters to ui ' , self.posx_antenna,self.ant_min_freq);
	self:update_create_model_settings_to_ui();
	self:update_antenna_settings_to_ui();
	self:update_solver_settings_to_ui();	
end
function RadomeSolver:update_parameters_from_ui()
	self:update_create_model_settings_from_ui();
	self:update_antenna_settings_from_ui();
	self:update_solver_settings_from_ui();
	--print( 'update parameters from ui ' , self.posx_antenna,self.ant_min_freq);
end
function RadomeSolver:update_create_model_settings_to_ui()
end
function RadomeSolver:update_create_model_settings_from_ui()
end
function RadomeSolver:update_antenna_settings_from_ui()
end
function RadomeSolver:update_antenna_settings_to_ui()
	print( 'invalid' )
end
function RadomeSolver:update_solver_settings_from_ui()
end
function RadomeSolver:update_solver_settings_to_ui()
end
function RadomeSolver:update_progress_bar()
end
--===========================================================================--

--===========================================================================--
--创建模型
function RadomeSolver:create_geometry(radome_path)
	--1.更新路径和材料设置
	print(radome_path)
	self.radome_path = radome_path;
	self:update_create_model_settings_from_ui()
	
	--2.参数检查
	assert( unity.file_exists( self.radome_path ), 'File do not exist： '..self.radome_path );
	assert( self.epsr > 0, 'epsr must grater than 0(介电常数需要大于0)： '..self.epsr );
	assert( self.tans >= 0, 'tans must grater than 0(损耗角正切需要大于或者等于0)： '..self.tans );
	if self.tans > 0 then 
		assert( self.tans_fit_freq > 0, 'fit frequecny must grater than 0： '..self.tans_fit_freq );
	end	
	
	--3.建立材料
	Geometry.clear_all();	
	local ufreq = Geometry.ref_units_freq();
	---------------------------
	--Materials：天线罩材料（外罩）
	if  self.tans > 1e-4 then
		self.material = Geometry.NewDebyeMaterial("dielectric")  --deby指的是色散模型
		self.material:fill_color(0.5,0.5,0)
		self.material:fit_tangent_delta_e(self.tans_fit_freq * Geometry.ref_units_freq(), self.epsr, self.tans)
	else
		self.material = Geometry.NewMaterial("dielectric")  --deby指的是色散模型
		self.material:fill_color(0.5,0.5,0)
		self.material:eps(self.epsr)
	end
	
	---------------------------

	--4.建立天线罩	
	--导入天线罩模型（外罩）：stp
	local w2 = Geometry.load_step({path = self.radome_path })
	w2:name( "天线罩" )
	w2:fill_color(0.6,0,0.8)
	w2:triangulate(1)
	w2:material( self.material )
	w2:soft_edge();
	self.loaded_shape = w2;
	Geometry.AddShape(self.loaded_shape)
	--------------------------- 
	
	Gallop.View:perspective_view();
end
--===========================================================================--



--===========================================================================--
--设置参数表
function RadomeSolver:create_paras()
	Geometry.Parameters.add_variable('scan_angle','0','天线扫描角（0°指向x轴负方向）');
	Geometry.Parameters.add_variable('posx_antenna','0','天线位置坐标');
	Geometry.Parameters.add_variable('radius_antenna','0','天线口径');
	Geometry.Parameters.add_variable('height','2');
	Geometry.Parameters.add_variable('dext','10','空气盒子与结构之间的间距');
	Geometry.Parameters.add_variable("polar_angle","0",'极化角,水平极化HH[0],垂直极化VV[90]')
end
function RadomeSolver:set_paras()
	
	local posx_antenna = self.posx_antenna;
	local radius_antenna = self.radius_antenna;
	assert(radius_antenna>1e-6,"天线口径设置错误");
	Geometry.Parameters.set_variable('posx_antenna',''..posx_antenna);
	Geometry.Parameters.set_variable('radius_antenna',''..radius_antenna);
	
	local lambda = c0/(self.mesh_freq * Geometry.ref_units_freq() ) / Geometry.ref_units_dim();
	Geometry.Parameters.set_variable('dext',''..lambda/2);

	local angs = self.scan_angles;
	local ang = angs[1];
	
	local polar_angles = nl.new_vector();
	if self.polarHH then 
		polar_angles:add( 0 );
	end
	if self.polarVV then 
		polar_angles:add( 90 );
	end
	local polarsStr = '';
	for kk,vv in pairs(polar_angles) do 
		polarsStr = polarsStr .. vv .. ';'
	end
	
	--if polar == 'HH' then
		--Geometry.Parameters.set_variable("polar_angle","90")
	--else
		--Geometry.Parameters.set_variable("polar_angle","0")
	--end
	Geometry.Parameters.set_variable('scan_angle',''.. self.min_angle);
	Geometry.Parameters.set_variable('polar_angle',''.. polar_angles[1]);
	
	
	Geometry.Parameters.change_variable_sweep_state("scan_angle",true)
	
	Geometry.Parameters.change_variable_sweep_values("scan_angle",'('.. self.min_angle .. ',' .. self.max_angle .. ',' .. self.step_angle .. ')')
		
	Geometry.Parameters.change_variable_sweep_state("polar_angle",true)
	
	Geometry.Parameters.change_variable_sweep_values("polar_angle",polarsStr)

	Geometry.Parameters.valuate()
	
	--print( Gallop.ParametersSweep.index );
	
	Gallop.ParametersSweep:change_solver_index('FDTD');
	
	Gallop.ParametersSweep:update_solver_list_to_ui();
	
	local preStr=[[	
-----------------------------------------------------------------------
--参数扫描前处理函数
function prerun(self)
	Gallop.RadomeSolverApp:initSolver();
end
-----------------------------------------------------------------------
]]

	local stepStr=[[	
-----------------------------------------------------------------------
--每一步扫描过程执行函数
function step(self,idx)
	--保存当前视图
	--self.view3D:saveSnapShot("view"..idx..".png")

	Gallop.RadomeSolverApp:preRunSolver(idx);

	--求解
	Gallop.Solver:runSolver();
	
	Gallop.RadomeSolverApp:postRunSolver(idx);	
		
end	
-----------------------------------------------------------------------
]]

	local postStr=[[	
-----------------------------------------------------------------------
--参数扫描后处理函数
function postrun(self)
end
-----------------------------------------------------------------------
]]

	Gallop.ParametersSweep:set_solver_settings( stepStr  );
	Gallop.ParametersSweep:set_pre_settings( preStr);
	Gallop.ParametersSweep:set_post_settings( postStr);
	Gallop.ParametersSweep:set_save_models(true);
	Gallop.ParametersSweep:set_export_gif(true);
	
	Geometry.synchronize();

end

--===========================================================================--

--===========================================================================--
function RadomeSolver:set_signal()
	self:update_antenna_settings_from_ui()
	local sig = Geometry.object( "pulse" , 'signal' );
	if sig == nil then sig = Geometry.NewGaussPulse("pulse") end
	sig:freq(self.ant_min_freq,self.ant_max_freq);
	self.sig = sig;
end
--===========================================================================--

--===========================================================================--
function RadomeSolver:set_nearfields()
	self:update_solver_settings_from_ui()
	if self.cal_nearfields then
		local probe = Geometry.object( "nearfields" , 'probe' );
		if probe == nil then 
			probe = Geometry.NewSurfaceProbe("nearfields") 
			Geometry.AddProbe(probe);
		end		
		local bbox = self.bbox;
		probe:freq( self.farbox_min_freq,self.farbox_max_freq,self.farbox_num_freq );
		probe:fields_types(ex,ey,ez)
		probe:pt1( bbox:xmin() .. '-dext', bbox:ymin() .. '-dext',   tostring((bbox:zmin()+bbox:zmax())/2) );
		probe:pt2( bbox:xmax() .. '+dext', bbox:ymax() .. '+dext',tostring((bbox:zmin()+bbox:zmax())/2) );
		if self.sig == nil then
			self:set_signal();
		end
		probe:signal(self.sig)
	else
		Geometry.remove_object( "nearfields" , 'probe' );
	end

end
--===========================================================================--

--===========================================================================--
function RadomeSolver:set_airbox()
	local bbox = Geometry.bbox()	
	assert( bbox:size() > 0, '没有任何几何结构' );
	
	local airbox = Geometry.object( "airbox" , 'geometry' );
	if airbox == nil then airbox = Geometry.NewBrick({name="airbox",is_fill=false}) end
	
	local objlist = Geometry.objects("shape");
	Geometry.remove_all_objects( "shapes" );
	Geometry.AddShape(airbox)
         for k,v in pairs ( objlist ) do
		if v:name() ~= "airbox"  then
			Geometry.AddShape(v)
		end
	end
	
	local objlist = Geometry.objects("shape");
	bbox = objlist[2]:bbox();
	for k=3,#objlist do 
		bbox:add( objlist[k]:bbox() );
	end
	airbox:pt1( bbox:xmin() ..'-dext', bbox:ymin()..'-dext', bbox:zmin() ..'-dext');
	airbox:pt2( bbox:xmax()..'+dext', bbox:ymax()..'+dext',bbox:zmax() ..'+dext');	
	self.bbox =  bbox;
end
--===========================================================================--

--===========================================================================--
function RadomeSolver:set_antenna()

	self:update_antenna_settings_from_ui()

	
	self:create_paras();
	self:set_paras();

	--print( "相扫？：",self.is_phasearray )
	
	Geometry.remove_all_objects('sources');
 	local srcbox = Geometry.create_cylinder();
	srcbox :radius('radius_antenna/2','0');
	srcbox :height('height');

	srcbox :rotate(yAxis,-90);
	if not self.is_phasearray then
		srcbox :rotate(zAxis,'-scan_angle');
	end
	srcbox :translate('posx_antenna','0','0')
	srcbox :is_fill(false);
	
	
	local plane = Geometry.create_plane()
	plane :rotate(yAxis,-90);
	if not self.is_phasearray then
		plane :rotate(zAxis,'-scan_angle');
	end
	plane :translate('posx_antenna','0','0')
	
	--local src = Geometry.GaussianBeamEx:new()
	self:set_signal();
	
	local src = Geometry.NewGaussianBeamEx({})	
	src :border(srcbox);
	src :signal(self.sig)
	src :direction(forward)
	src :width_u('radius_antenna/2*0.9')
	src :width_v('radius_antenna/2*0.9')
	src :rotate(zAxis,-90);
	src :rotate(yAxis,-90);
	src :rotate(zAxis,'-scan_angle');
	src :translate('posx_antenna','0','0')
	
	src :point_size(2)
	src :add_cut_plane(plane);
	src :is_point(true)
	src :outline_width(1.5)
	src :polarization_angle('polar_angle')
	src :theta_ref_freq((self.ant_min_freq+self.ant_max_freq)/2.0)
	Geometry.AddSrc(src)		
	self.bbox = src:bbox()
end
--===========================================================================--

--===========================================================================--
--设置抓屏器
function RadomeSolver:set_grabber()
	local bbox = Geometry.bbox()	
	assert( bbox:size() > 0, '没有任何几何结构' );
	
	local grabber = Geometry.object( "grabber" , 'grabber' );
	if grabber == nil then 
		grabber = Geometry.NewGrabber({name="grabber"}) 
		cm = Geometry.CGLColorMap:new("Jet");
		cm:min(0)
		
		grabber:add_slice(cm,scalar_e,xy,bbox:zcenter(),0,false,'xy 平面')
		grabber:add_slice(cm,scalar_e,zx,bbox:ycenter(),0,false,'xz 平面')
		Geometry.AddTransientFieldsGrabber(grabber)	
	end
end
--===========================================================================--

--===========================================================================--
function RadomeSolver:set_solver()
	self:update_solver_settings_from_ui()
	
	self:create_paras();
	--self:set_paras();
	
	Geometry.Solver.threads( self.num_threads );
	Geometry.Solver.accuracy( self.accuracy );
	Geometry.Solver.max_steps( self.max_steps );
	
	assert( self.farbox_min_freq <= self.farbox_max_freq, '求解频率设置错误')
	assert( self.mesh_freq >= 1e-6, '剖分频率设置错误')
	Geometry.Mesh.mesh_freq(self.mesh_freq)
	local lambda = c0/(self.mesh_freq * Geometry.ref_units_freq() ) / Geometry.ref_units_dim();
	Geometry.remove_all_objects('farbox');
	Geometry.Parameters.set_variable('dext',''..lambda/2);
	
	self:set_airbox();
	
	Geometry.Boundaries.boundary(UPML,UPML,UPML,UPML,UPML,UPML);
	Geometry.ref_units(mm,GHz,ns)
	local bbox = self.bbox;
	
	assert( bbox:size() > 0, '没有任何几何结构' );

	local farbox = Geometry.object( "fbox" );
	if farbox == nil then 
		farbox = Geometry.NewNear2FarFieldBox("fbox") 
		Geometry.AddProbe(farbox);
	end
	farbox:freq( self.farbox_min_freq,self.farbox_max_freq,self.farbox_num_freq );
	farbox:pt1( bbox:xmin() .. '-dext', bbox:ymin() .. '-dext', bbox:zmin() .. '-dext' );
	farbox:pt2( bbox:xmax() .. '+dext', bbox:ymax() .. '+dext',bbox:zmax() .. '+dext' );
	
	farbox:angle_theta( 90,90,1);
	farbox:angle_phi( '180-scan_angle-30','180-scan_angle+30',301);


	if self.sig == nil then
		self:set_signal();
	end
	farbox:signal(self.sig)
	
	
	self:set_nearfields();
	
	self:set_grabber()
	
end
--===========================================================================--

--===========================================================================--
function RadomeSolver:read_from_files(str)
	print(str)
	x = nl.LinearSpace(0,1,10);
	y = nl.LinearSpace(0,1,10);
	self.tranView:addGraph(x,y);
	self.tranView:rescaleAxes();
end
--===========================================================================--

--===========================================================================--
function RadomeSolver:run_pond()
	print('run pond')
	--assert(self.binit, '模型未初始化' )		 
	Gallop.ParametersSweep:run_sweep()
end
--===========================================================================--

function RadomeSolver:initSolver()
	if self.bseView then
		self.bseView:clearGraphs()
		self.bseView:rescaleAxes();
	end
	
	if self.bssView then
		self.bssView:clearGraphs()
		self.bssView:rescaleAxes()
	end
	
	if self.tranView then
		self.tranView:clearGraphs()
		self.tranView:rescaleAxes()
	end
	
	if self.farfieldView then
		self.farfieldView:clearGraphs()
		self.farfieldView:rescaleAxes()
	end
	self.results={};	
end

--===========================================================================--
function RadomeSolver:preRunSolver(idx)
	print( "============Start============" )
	print( "扫参:"..idx )
	function index( names, name )
		for k,v in pairs(names) do
			if v==name then
				return k
			end
		end
	end
	
	local names = Gallop.ParametersSweep.names;
	local polar_angle = 0;
	if index(names,'polar_angle') ~= nil then
		polar_angle = Gallop.ParametersSweep.paraList[idx][index(names,'polar_angle')];
	end
	
	local polar = 'VV'
	if nl.abs( polar_angle - 90) < 1e-3 then
		polar = 'HH';
	end
	local scan_angle = Gallop.ParametersSweep.paraList[idx][index(names,'scan_angle')];
        self.curr_lable = 'angle-'..scan_angle..' polar-'..polar;
	
	local grabber = Geometry.objects('Grabber');
	if grabber ~= nil and #grabber > 0 then
		grabber[1]:name('Grabber('.. self.curr_lable .. ')');
	end
end
--===========================================================================--

--===========================================================================--
function RadomeSolver:postRunSolver(idx)
	print( "post calculate finished, idx = [".. idx .. "]\n" )	
	self:update_progress_bar( idx / #Gallop.ParametersSweep.paraList * 100 );
	print( "==============End=============" )
end
--===========================================================================--



--===========================================================================--
--读路径文件
function RadomeSolver:read_background_from_files(path)
	print( path )
	local data = self.data;
	self:read_from_files(path);
	self.background = self.data;
	self.data = data;
end
--===========================================================================--

--===========================================================================--
--读路径文件
function RadomeSolver:read_from_files(path)
	--print( path )
	--第一步：解析路径文件名称
	self:parser_files(path);
	
	--第二步：重新索引
	self:reindex()
	
	--第三步：获取文件包含的角度和极化
	self:cal_info()
	
	--第四步：读取文件
	for k,v in pairs( self.results_files ) do 
		local phi, etheta_trans,ephi_trans = self:read_file(v.path)
		v.phi = phi;
		v.etheta_trans = etheta_trans;		
		v.ephi_trans = ephi_trans;		
	end
	
	self.data = nl.new_eigen_tensor3( #self.polars, #self.angles, self.nfreq )
	for kp,vp in pairs( self.polars) do 
		for ka,va in pairs( self.angles) do 
			if vp == '0.0' then
				for kf,vf in pairs( self.results_files[ vp..va ].ephi_trans ) do 
					self.data[{kp,ka,kf}] = vf;
				end
			end
			if vp == '90.0' then
				for kf,vf in pairs( self.results_files[ vp..va ].etheta_trans ) do 
					self.data[{kp,ka,kf}] = vf;
				end
			end
		end
	end
	--table.sort( self.results_files, function(a,b) return a.angle<b.angle end );
	self:plot_trans();
end
--===========================================================================--

--===========================================================================--
--读文件角度和极化信息
function RadomeSolver:cal_info()
	local get_key = function( data, name )
		local keys = {};
		for k,v in pairs( data ) do 
			keys[v[name]] = true;		
		end
		local res ={};
		for k,v in pairs(keys) do 
			table.insert(res,k);
		end
		return res;
	end
	self.polars = get_key(self.results_files,'polar');
	self.angles = get_key(self.results_files,'angle');
	table.sort( self.polars, function(a,b) return tonumber(a)<tonumber(b) end );
	table.sort( self.angles, function(a,b) return tonumber(a)<tonumber(b) end );

end
--===========================================================================--

--===========================================================================--
--重新索引
function RadomeSolver:reindex()
	local results_files = {};
	for k,v in pairs( self.results_files ) do 
		results_files[ v.key ] = v;		
	end
	 self.results_files = results_files;
end
--===========================================================================--

--===========================================================================--
--读路径文件
function RadomeSolver:plot_trans()
	self.tranView:clearGraphs();
	local angles = nl.new_vector( #self.angles );
	for ka,va in pairs( self.angles) do 
		angles[ka] = tonumber(va);
	end

	for kp,vp in pairs( self.polars) do 
		for kf=1,self.nfreq do 
			local trans = nl.new_vector( #self.angles );
			for ka,va in pairs( self.angles) do 
				trans[ka] = self.data[{kp,ka,kf}];
				if self.background then 
					trans[ka] = self.data[{kp,ka,kf}]/self.background[{kp,ka,kf}];
					trans[ka]  = trans[ka] * trans[ka]  * 100;
				end
			end
			self.tranView:addGraph(angles,trans,'polar angle:'..vp..' Freq:'..kf)
		end
	end
end
--===========================================================================--

--===========================================================================--
--解析文件路径：获得符合命名规则的文件名称
function RadomeSolver:parser_files(path)
	local strScanAngle = "scan_angle%((.-)%)";
	local strPolarAngle = "polar_angle%((.-)%)";
	local strName = "%[name%-%((.-)%)%]";
	self.results_files={}
	for file in lfs.dir(path) do
	        if file ~= "." and file ~= ".." then
			local attr = lfs.attributes(path.."/"..file);
			if attr.mode ~=  "directory" then
				local angle,polar,name = string.match(file, strScanAngle .. '.-' ..strPolarAngle .. '.-' .. strName ..'.-%.txt')
				if polar and angle and name=='fbox' then
					--print( polar,angle ,name)
					table.insert( self.results_files, { key = polar..angle, 
											path=path.."/"..file, 
											polar=polar;angle=angle; 
											numangle  = tonumber(angle);
											name = name} )
				end
			end
		end
	end
end
--===========================================================================--

--===========================================================================--
--读取文件
function RadomeSolver:read_file(path)
	local colvector = function( data, ncol )
		local dvec = nl.new_vector( #data );
		for kk=1,#dvec do 
			dvec[kk] = data[kk][ncol];
		end
		return dvec;
	end
	local data = unity.file_to_num_table(path,'\t',8,1);
	--print( #data,#data[1] )
	local nrow = #data;
	local nfreq = math.floor((#data[1] - 3)/8);
	
	local phi = colvector( data, 2);
	local etheta_trans = {};
	local ephi_trans = {};
	for kk=1,nfreq do
		local farfields = colvector( data, 2+(kk-1)*8+1 );
		etheta_trans[kk] = nl.max( farfields );
		local farfields = colvector( data, 2+(kk-1)*8+3 );
		ephi_trans[kk] = nl.max( farfields );
	end
	self.nfreq = nfreq;
	return phi, etheta_trans,ephi_trans;	
	--print_table( data[1],#data[1] )
end
--===========================================================================--

return RadomeSolver;
